Взять дневные или недельные цены акций 5-7 компаний за период 2015-2021 годы включительно. Из них 2 компании должны быть те, кто имеет привилегированные акции.
Построить одномерные модели временных рядов, оценить волатильность.
Построить модель векторной авторегрессии для этих акций, проверить каждое уравнение VAR модели на структурный сдвиг в период начала пандемии (в США - 24 февраля) и в момент выхода новости о создании вакцины (в США - 4 ноября).
Выяснить, наблюдается ли коинтеграция между рассматриваемыми акциями? Продемонстрировать, как коинтеграция влияет на прогноз.
Для расчета ожидаемой доходности финансовых инструментов использовать модель CAPM.
Для данного исследования мы выбрали акции следующих компаний:
ПАО “ТАТНЕФТЬ” (TATN_p) - привилегированные акции
ПАО “СУРГУТНЕФТЕГАЗ” (SNGS_p) - привилегированные акции
ОАО АК “АЛРОСА” (ALRS)
ПАО “ГАЗПРОМ НЕФТЬ” (SIBN)
ПАО “НОВОЛИПЕЦКИЙ МЕТАЛЛУРГИЧЕСКИЙ КОМБИНАТ” (NLMK)
ПАО “КАМАЗ” (KMAZ)
ПАО “АЭРОФЛОТ-РОССИЙСКИЕ АВИАЛИНИИ” (AFLT)
Мы взяли дневные данные с 05.01.2015 по 30.12.2021 (исключая выходные на бирже). Так они выглядят в исходнои виде:
head(data, 5)
tail(data, 5)
time = data$date
time = as.Date(time, "%m/%d/%y")
plot(time, data$TATN_p)
plot(time, data$SNGS_p)
plot(time, data$ALRS)
plot(time, data$SIBN)
plot(time, data$NLMK)
plot(time, data$KMAZ)
plot(time, data$AFLT)
На графиках сезонность не наблюдается, что ожидаемо, так как цены на акции меняются стихийно, вследствие прочих внешних факторов.
Критическое значение для проверки на 1% уровне значимости = -3.43
#Pacf(data$TATN_p) #lags = 1
stTANT = ur.df(data$TATN_p, type="drift", lags = 1, selectlags = "Fixed")
stTANT@teststat
## tau2 phi1
## statistic -1.556173 1.577026
Значение t статистики -1.556173 > -3.43.
#Pacf(data$SNGS_p) #lags = 1
stSNGS = ur.df(data$SNGS_p, type="drift", lags = 1, selectlags = "Fixed")
stSNGS@teststat
## tau2 phi1
## statistic -2.913332 4.291275
Значение t статистики -2.913332 > -3.43.
#Pacf(data$ALRS) #lags = 1
stALRS = ur.df(data$ALRS, type="drift", lags = 1, selectlags = "Fixed")
stALRS@teststat
## tau2 phi1
## statistic -1.520594 1.532856
Значение t статистики -1.520594 > -3.43.
#Pacf(data$SIBN) #lags = 1
stSIBN = ur.df(data$SIBN, type="drift", lags = 1, selectlags = "Fixed")
stSIBN@teststat
## tau2 phi1
## statistic -0.001731526 1.700032
Значение t статистики -0.001731526 > -3.43.
#Pacf(data$NLMK) #lags = 1
stNLMK = ur.df(data$NLMK, type="drift", lags = 1, selectlags = "Fixed")
stNLMK@teststat
## tau2 phi1
## statistic -0.9758007 1.515227
Значение t статистики -0.9758007 > -3.43.
#Pacf(data$KMAZ) #lags = 1
stKMAZ = ur.df(data$KMAZ, type="drift", lags = 1, selectlags = "Fixed")
stKMAZ@teststat
## tau2 phi1
## statistic -0.9685645 1.289144
Значение t статистики -0.9685645 > -3.43.
#Pacf(data$AFLT) #lags = 1
stAFLT = ur.df(data$AFLT, type="drift", lags = 1, selectlags = "Fixed")
stAFLT@teststat
## tau2 phi1
## statistic -1.397104 1.014872
Значение t статистики -1.397104 > -3.43.
Гипотеза о нестационарности рядов не отклоняется => ряды нестационарны, применяем тест Энгла-Гренджера на коинтеграцию для начальных данных, а для построения ARMA-модели будем использовать разности.
d1TATN = diff(data$TATN_p, differences=1)
d1SNGS = diff(data$SNGS_p, differences=1)
d1ALRS = diff(data$ALRS, differences=1)
d1SIBN = diff(data$SIBN, differences=1)
d1NLMK = diff(data$NLMK, differences=1)
d1KMAZ = diff(data$KMAZ, differences=1)
d1AFLT = diff(data$AFLT, differences=1)
#Pacf(d1TATN) #lags = 2
std1TATN = ur.df(d1TATN, type="drift", lags = 2, selectlags = "Fixed")
std1TATN@teststat
## tau2 phi1
## statistic -21.80503 237.73
Значение t статистики -21.80503 < -3.43.
#Pacf(d1SNGS) #lags = 7
std1SNGS = ur.df(d1SNGS, type="drift", lags = 7, selectlags = "Fixed")
std1SNGS@teststat
## tau2 phi1
## statistic -14.47239 104.7288
Значение t статистики -14.47239 < -3.43.
#Pacf(d1ALRS) #lags = 27
std1ALRS = ur.df(d1ALRS, type="drift", lags = 27, selectlags = "Fixed")
std1ALRS@teststat
## tau2 phi1
## statistic -7.792801 30.36435
Значение t статистики -7.792801 < -3.43.
#Pacf(d1SIBN) #lags = 1
std1SIBN = ur.df(d1SIBN, type="drift", lags = 1, selectlags = "Fixed")
std1SIBN@teststat
## tau2 phi1
## statistic -28.52974 406.9731
Значение t статистики -28.52974 < -3.43.
#Pacf(d1NLMK) #lags = 7
std1NLMK = ur.df(d1NLMK, type="drift", lags = 7, selectlags = "Fixed")
std1NLMK@teststat
## tau2 phi1
## statistic -13.70413 93.90189
Значение t статистики -13.70413 < -3.43.
#Pacf(d1KMAZ) #lags = 2
std1KMAZ = ur.df(d1KMAZ, type="drift", lags = 2, selectlags = "Fixed")
std1KMAZ@teststat
## tau2 phi1
## statistic -26.38805 348.1647
Значение t статистики -26.38805 < -3.43.
#Pacf(d1AFLT) #lags = 1
std1AFLT = ur.df(d1AFLT, type="drift", lags = 1, selectlags = "Fixed")
std1AFLT@teststat
## tau2 phi1
## statistic -28.56271 407.9146
Значение t статистики -28.56271 < -3.43.
Гипотеза о нестационарности отклоняется на 1% уровне => разности стационарны, можно строить ARMA модель на первых разностях.
lag = ln(1765) = 7
Параметры ARMA(p, k) выбирались на основе сравнения BIC и AIC нескольких моделей ARMA.
#eacf(d1TATN) # (p, k) = (0,2)
ARMAmodelTANT = Arima(data$TATN_p, c(0,1,2), include.constant=TRUE, method = c("CSS-ML"))
#Acf(residuals(ARMAmodelTANT)) # автокорреляция на графике
Box.test(residuals(ARMAmodelTANT), lag = 7, type = c("Ljung-Box"), fitdf = 2) # тест на автокорреляцию
##
## Box-Ljung test
##
## data: residuals(ARMAmodelTANT)
## X-squared = 7.0018, df = 5, p-value = 0.2205
shapiro.test(residuals(ARMAmodelTANT)) # проверка на нормальное распределение
##
## Shapiro-Wilk normality test
##
## data: residuals(ARMAmodelTANT)
## W = 0.83457, p-value < 2.2e-16
ARMAmodelSNGS = Arima(data$SNGS_p, c(0,1,1), include.constant=TRUE, method = c("CSS-ML"))
Acf(residuals(ARMAmodelSNGS))
Box.test(residuals(ARMAmodelSNGS), lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(ARMAmodelSNGS)
## X-squared = 9.002, df = 6, p-value = 0.1735
shapiro.test(residuals(ARMAmodelSNGS))
##
## Shapiro-Wilk normality test
##
## data: residuals(ARMAmodelSNGS)
## W = 0.81328, p-value < 2.2e-16
ARMAmodelALRS = Arima(data$ALRS, c(0,1,1), include.constant=TRUE, method = c("CSS-ML"))
Acf(residuals(ARMAmodelALRS))
Box.test(residuals(ARMAmodelALRS), lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(ARMAmodelALRS)
## X-squared = 3.8173, df = 6, p-value = 0.7014
shapiro.test(residuals(ARMAmodelALRS))
##
## Shapiro-Wilk normality test
##
## data: residuals(ARMAmodelALRS)
## W = 0.97255, p-value < 2.2e-16
ARMAmodelSIBN = Arima(data$SIBN, c(0,1,4), include.constant=TRUE, method = c("CSS-ML"))
Acf(residuals(ARMAmodelSIBN))
Box.test(residuals(ARMAmodelSIBN), lag = 7, type = c("Ljung-Box"), fitdf = 4)
##
## Box-Ljung test
##
## data: residuals(ARMAmodelSIBN)
## X-squared = 2.6499, df = 3, p-value = 0.4488
shapiro.test(residuals(ARMAmodelSIBN))
##
## Shapiro-Wilk normality test
##
## data: residuals(ARMAmodelSIBN)
## W = 0.86721, p-value < 2.2e-16
ARMAmodelNLMK = Arima(data$NLMK, c(0,1,1), include.constant=TRUE, method = c("CSS-ML"))
Acf(residuals(ARMAmodelNLMK))
Box.test(residuals(ARMAmodelNLMK), lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(ARMAmodelNLMK)
## X-squared = 6.3029, df = 6, p-value = 0.3901
shapiro.test(residuals(ARMAmodelNLMK))
##
## Shapiro-Wilk normality test
##
## data: residuals(ARMAmodelNLMK)
## W = 0.96844, p-value < 2.2e-16
ARMAmodelKMAZ = Arima(data$KMAZ, c(0,1,4), include.constant=TRUE, method = c("CSS-ML"))
Acf(residuals(ARMAmodelKMAZ))
Box.test(residuals(ARMAmodelKMAZ), lag = 7, type = c("Ljung-Box"), fitdf = 4)
##
## Box-Ljung test
##
## data: residuals(ARMAmodelKMAZ)
## X-squared = 2.3316, df = 3, p-value = 0.5065
shapiro.test(residuals(ARMAmodelKMAZ))
##
## Shapiro-Wilk normality test
##
## data: residuals(ARMAmodelKMAZ)
## W = 0.6537, p-value < 2.2e-16
ARMAmodelAFLT = Arima(data$AFLT, c(0,1,1), include.constant=TRUE, method = c("CSS-ML"))
Acf(residuals(ARMAmodelAFLT))
Box.test(residuals(ARMAmodelAFLT), lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(ARMAmodelAFLT)
## X-squared = 6.212, df = 6, p-value = 0.3999
shapiro.test(residuals(ARMAmodelAFLT))
##
## Shapiro-Wilk normality test
##
## data: residuals(ARMAmodelAFLT)
## W = 0.89919, p-value < 2.2e-16
Box-Ljung test: p-value > 0.05 - автокорреляция остатков отсутствует на уровне значимости 5%, на графиках можно увидеть аналогичный результат.
Shapiro-Wilk normality test: p-value < 2.2e-16 - гипотеза о нормальном распределении отвергается.
Параметры и тип модели GARCH выбирались опытным путем через сравнение моделей по показателям LogLikelihood и Information Criteria.
При расчетах используется первая разность.
Также каждая модель проверяется на автокорреляцию стандатизированных остатков и остатков в квадрате с помощью Box-Ljung теста.
spec_TATN = ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 2), include.mean = TRUE), distribution.model = "norm")
garch.TATN = ugarchfit(spec_TATN, d1TATN)
Box.test(residuals(garch.TATN, standardize="TRUE"), lag = 7, type = c("Ljung-Box"), fitdf = 2)
##
## Box-Ljung test
##
## data: residuals(garch.TATN, standardize = "TRUE")
## X-squared = 4.9887, df = 5, p-value = 0.4173
Box.test(residuals(garch.TATN, standardize="TRUE")^2, lag = 7, type = c("Ljung-Box"), fitdf = 2)
##
## Box-Ljung test
##
## data: residuals(garch.TATN, standardize = "TRUE")^2
## X-squared = 5.668, df = 5, p-value = 0.3399
spec_SNGS = ugarchspec(variance.model = list(model = 'sGARCH',garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 1), include.mean = TRUE), distribution.model = "std")
garch.SNGS = ugarchfit(spec_SNGS, d1SNGS)
Box.test(residuals(garch.SNGS, standardize="TRUE"), lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(garch.SNGS, standardize = "TRUE")
## X-squared = 3.561, df = 6, p-value = 0.7358
Box.test(residuals(garch.SNGS, standardize="TRUE")^2, lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(garch.SNGS, standardize = "TRUE")^2
## X-squared = 0.84517, df = 6, p-value = 0.9908
spec_ALRS = ugarchspec(variance.model = list(model = 'sGARCH',garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 1), include.mean = TRUE), distribution.model = "norm")
garch.ALRS = ugarchfit(spec_ALRS, d1ALRS)
Box.test(residuals(garch.ALRS, standardize="TRUE"), lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(garch.ALRS, standardize = "TRUE")
## X-squared = 6.9212, df = 6, p-value = 0.3282
Box.test(residuals(garch.ALRS, standardize="TRUE")^2, lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(garch.ALRS, standardize = "TRUE")^2
## X-squared = 7.4364, df = 6, p-value = 0.2824
spec_SIBN = ugarchspec(variance.model = list(model = 'sGARCH',garchOrder = c(1, 3)),
mean.model = list(armaOrder = c(0, 4), include.mean = TRUE), distribution.model = "std")
garch.SIBN = ugarchfit(spec_SIBN, d1SIBN)
Box.test(residuals(garch.SIBN, standardize="TRUE"), lag = 7, type = c("Ljung-Box"), fitdf = 4)
##
## Box-Ljung test
##
## data: residuals(garch.SIBN, standardize = "TRUE")
## X-squared = 5.5013, df = 3, p-value = 0.1386
Box.test(residuals(garch.SIBN, standardize="TRUE")^2, lag = 7, type = c("Ljung-Box"), fitdf = 4)
##
## Box-Ljung test
##
## data: residuals(garch.SIBN, standardize = "TRUE")^2
## X-squared = 4.5245, df = 3, p-value = 0.2101
spec_NLMK = ugarchspec(variance.model = list(model = 'sGARCH',garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 1), include.mean = TRUE), distribution.model = "std")
garch.NLMK = ugarchfit(spec_NLMK, d1NLMK)
Box.test(residuals(garch.NLMK, standardize="TRUE"), lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(garch.NLMK, standardize = "TRUE")
## X-squared = 3.7148, df = 6, p-value = 0.7152
Box.test(residuals(garch.NLMK, standardize="TRUE")^2, lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(garch.NLMK, standardize = "TRUE")^2
## X-squared = 6.6174, df = 6, p-value = 0.3577
spec_KMAZ = ugarchspec(variance.model = list(model = 'sGARCH',garchOrder = c(1, 3)),
mean.model = list(armaOrder = c(0, 4), include.mean = TRUE), distribution.model = "std")
garch.KMAZ = ugarchfit(spec_KMAZ, d1KMAZ)
Box.test(residuals(garch.KMAZ, standardize="TRUE"), lag = 7, type = c("Ljung-Box"), fitdf = 4)
##
## Box-Ljung test
##
## data: residuals(garch.KMAZ, standardize = "TRUE")
## X-squared = 13.125, df = 3, p-value = 0.004373
Box.test(residuals(garch.KMAZ, standardize="TRUE")^2, lag = 7, type = c("Ljung-Box"), fitdf = 4)
##
## Box-Ljung test
##
## data: residuals(garch.KMAZ, standardize = "TRUE")^2
## X-squared = 0.52356, df = 3, p-value = 0.9137
spec_AFLT = ugarchspec(variance.model = list(model = 'sGARCH',garchOrder = c(1, 3)),
mean.model = list(armaOrder = c(0, 1), include.mean = TRUE), distribution.model = "std")
garch.AFLT = ugarchfit(spec_AFLT, d1AFLT)
Box.test(residuals(garch.AFLT, standardize="TRUE"), lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(garch.AFLT, standardize = "TRUE")
## X-squared = 11.503, df = 6, p-value = 0.07403
Box.test(residuals(garch.AFLT, standardize="TRUE")^2, lag = 7, type = c("Ljung-Box"), fitdf = 1)
##
## Box-Ljung test
##
## data: residuals(garch.AFLT, standardize = "TRUE")^2
## X-squared = 2.0451, df = 6, p-value = 0.9155
Во всех случаях в стандатизированных остатках в квадрате отсутсвует автокорреляция, что является признаком хорошей модели.
На коинтеграцию в тесте Энгла-Гренджера проверяются все нестационарные данные. В нашем случае - все начальные ряды.
Критическое значение на уровне 5% без тренда = -3.34, на уровне 1% без тренда = -3.90. Данные коинтегриваны при значении t-статистики < критического.
coint_TATN_SNGS = lm(data$TATN_p ~ data$SNGS_p)
coint_TATN_ALRS = lm(data$TATN_p ~ data$ALRS)
coint_TATN_SIBN = lm(data$TATN_p ~ data$SIBN)
coint_TATN_NLMK = lm(data$TATN_p ~ data$NLMK)
coint_TATN_KMAZ = lm(data$TATN_p ~ data$KMAZ)
coint_TATN_AFLT = lm(data$TATN_p ~ data$AFLT)
ur.df(coint_TATN_SNGS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.5013 1.4453
ur.df(coint_TATN_ALRS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.446 1.1099
ur.df(coint_TATN_SIBN$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: 0.2291 0.4253
ur.df(coint_TATN_NLMK$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.4745 1.0874
ur.df(coint_TATN_KMAZ$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.3097 0.9066
ur.df(coint_TATN_AFLT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.4967 1.4806
С TATN_p коинтеграции нет.
coint_SNGS_TANT = lm(data$SNGS_p ~ data$TATN_p)
coint_SNGS_ALRS = lm(data$SNGS_p ~ data$ALRS)
coint_SNGS_SIBN = lm(data$SNGS_p ~ data$SIBN)
coint_SNGS_NLMK = lm(data$SNGS_p ~ data$NLMK)
coint_SNGS_KMAZ = lm(data$SNGS_p ~ data$KMAZ)
coint_SNGS_AFLT = lm(data$SNGS_p ~ data$AFLT)
ur.df(coint_SNGS_TANT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.8917 4.2211
ur.df(coint_SNGS_ALRS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.895 4.23
ur.df(coint_SNGS_SIBN$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.8964 4.2235
ur.df(coint_SNGS_NLMK$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.897 4.2225
ur.df(coint_SNGS_KMAZ$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.9049 4.2601
ur.df(coint_SNGS_AFLT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -5.4669 15.0145
С SNGS_p коинтегрированна только переменная AFLT.
coint_ALRS_TANT = lm(data$ALRS ~ data$TATN_p)
coint_ALRS_SNGS = lm(data$ALRS ~ data$SNGS_p)
coint_ALRS_SIBN = lm(data$ALRS ~ data$SIBN)
coint_ALRS_NLMK = lm(data$ALRS ~ data$NLMK)
coint_ALRS_KMAZ = lm(data$ALRS ~ data$KMAZ)
coint_ALRS_AFLT = lm(data$ALRS ~ data$AFLT)
ur.df(coint_ALRS_TANT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.4831 1.3459
ur.df(coint_ALRS_SNGS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.4893 1.4717
ur.df(coint_ALRS_SIBN$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.2661 2.59
ur.df(coint_ALRS_NLMK$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.7277 3.7439
ur.df(coint_ALRS_KMAZ$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -3.0689 4.7165
ur.df(coint_ALRS_AFLT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.3729 1.3105
С ALRS коинтеграции нет.
coint_SIBN_TANT = lm(data$SIBN ~ data$TATN_p)
coint_SIBN_SNGS = lm(data$SIBN ~ data$SNGS_p)
coint_SIBN_ALRS = lm(data$SIBN ~ data$ALRS)
coint_SIBN_NLMK = lm(data$SIBN ~ data$NLMK)
coint_SIBN_KMAZ = lm(data$SIBN ~ data$KMAZ)
coint_SIBN_AFLT = lm(data$SIBN ~ data$AFLT)
ur.df(coint_SIBN_TANT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: 1.1736 1.5028
ur.df(coint_SIBN_SNGS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: 0.0343 1.5387
ur.df(coint_SIBN_ALRS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.4087 1.2166
ur.df(coint_SIBN_NLMK$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.3987 1.157
ur.df(coint_SIBN_KMAZ$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.9036 4.2163
ur.df(coint_SIBN_AFLT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -0.0155 1.6982
С SIBN коинтеграции нет.
coint_NLMK_TANT = lm(data$NLMK ~ data$TATN_p)
coint_NLMK_SNGS = lm(data$NLMK ~ data$SNGS_p)
coint_NLMK_ALRS = lm(data$NLMK ~ data$ALRS)
coint_NLMK_SIBN = lm(data$NLMK ~ data$SIBN)
coint_NLMK_KMAZ = lm(data$NLMK ~ data$KMAZ)
coint_NLMK_AFLT = lm(data$NLMK ~ data$AFLT)
ur.df(coint_NLMK_TANT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.0397 0.8438
ur.df(coint_NLMK_SNGS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -0.9367 1.2336
ur.df(coint_NLMK_ALRS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.471 3.0671
ur.df(coint_NLMK_SIBN$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.9031 1.8112
ur.df(coint_NLMK_KMAZ$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.4753 3.0933
ur.df(coint_NLMK_AFLT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -0.9753 1.5147
С NLMK коинтеграции нет.
coint_KMAZ_TANT = lm(data$KMAZ ~ data$TATN_p)
coint_KMAZ_SNGS = lm(data$KMAZ ~ data$SNGS_p)
coint_KMAZ_ALRS = lm(data$KMAZ ~ data$ALRS)
coint_KMAZ_SIBN = lm(data$KMAZ ~ data$SIBN)
coint_KMAZ_NLMK = lm(data$KMAZ ~ data$NLMK)
coint_KMAZ_AFLT = lm(data$KMAZ ~ data$AFLT)
ur.df(coint_KMAZ_TANT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.0175 1.0046
ur.df(coint_KMAZ_SNGS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -0.9542 1.2625
ur.df(coint_KMAZ_ALRS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.7516 3.9455
ur.df(coint_KMAZ_SIBN$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -3.1226 4.9635
ur.df(coint_KMAZ_NLMK$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -2.416 3.1114
ur.df(coint_KMAZ_AFLT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -0.8838 1.2057
С KMAZ коинтеграции нет.
coint_AFLT_TANT = lm(data$AFLT ~ data$TATN_p)
coint_AFLT_SNGS = lm(data$AFLT ~ data$SNGS_p)
coint_AFLT_ALRS = lm(data$AFLT ~ data$ALRS)
coint_AFLT_SIBN = lm(data$AFLT ~ data$SIBN)
coint_AFLT_NLMK = lm(data$AFLT ~ data$NLMK)
coint_AFLT_KMAZ = lm(data$AFLT ~ data$KMAZ)
ur.df(coint_AFLT_TANT$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.3352 0.9198
ur.df(coint_AFLT_SNGS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -4.8218 11.7055
ur.df(coint_AFLT_ALRS$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.1566 0.6728
ur.df(coint_AFLT_SIBN$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.4083 1.0331
ur.df(coint_AFLT_NLMK$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.3966 1.0141
ur.df(coint_AFLT_KMAZ$residuals, type="drift", lags = 1, selectlags = "Fixed")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -1.2597 0.8071
С AFLT коинтегрированна только переменная SNGS_p.
Из всех данных коинтегрированными на 1% уровне является только пара SNGS_p и AFLT. Данную коинтеграцию можно использовать в качестве основной модели для построения совместного прогноза.
Для остальных переменных - переходить на первую разность и строить вектор авторегрессионной модели.
Проверим все акции на коинтеграцию:
df_all = data %>% dplyr::select(-date)
coint_jo_all = ca.jo(df_all, ecdet = "const", type = "eigen", K=2, spec = "transitory", season = NULL)
jo = data.frame(coint_jo_all@teststat, coint_jo_all@cval)
jo
На уровне r = 0: 43.147034 < 43.25 => нулевая гипотеза (коинтеграция между 7 переменными отсутствует) не отвергается, то есть коинтеграции между выбранными нами акциями не неблюдается.
Таким образом, для совместного прогноза берем разности по всем данным и строим VAR model.
Целевая переменная - TATN_p.
ccf(d1TATN, d1SNGS, lag.max = 20, type = c("correlation"), plot = TRUE)
grangertest(d1SNGS, d1TATN, order = 1)
p-value = 0.4858 => влияния нет
ccf(d1TATN, d1ALRS, lag.max = 35, type = c("correlation"), plot = TRUE)
grangertest(d1ALRS, d1TATN, order = 31)
p-value = 0.3605 => влияния нет
ccf(d1TATN, d1SIBN, lag.max = 30, type = c("correlation"), plot = TRUE)
grangertest(d1SIBN, d1TATN, order = 29)
p-value = 2.718e-06 => влияние есть
ccf(d1TATN, d1NLMK, lag.max = 25, type = c("correlation"), plot = TRUE)
grangertest(d1NLMK, d1TATN, order = 19)
p-value = 0.009423 => влияние есть
ccf(d1TATN, d1KMAZ, lag.max = 30, type = c("correlation"), plot = TRUE)
grangertest(d1KMAZ, d1TATN, order = 25)
p-value = 0.005577 => влияние есть
ccf(d1TATN, d1AFLT, lag.max = 30, type = c("correlation"), plot = TRUE)
grangertest(d1AFLT, d1TATN, order = 22)
p-value = 0.05897 => влияние есть
Наиболее сильное влияние имеют SIBN, KMAZ, NLMK и AFLT.
Целевая переменная анализа - привилегированные акции компании ПАО “Татнефть”.
df = data.frame(d1TATN, d1SIBN)
VARselect(df, lag.max = 23, type = 'const')$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 10 2 1 10
var = VAR(df, p = 10, type = 'const')
Hosking(var, lags=1.5*var$p)
## lags statistic df p-value
## 15 20.99677 20 0.3973226
LiMcLeod(var, lags=1.5*var$p)
## lags statistic df p-value
## 15 21.11511 20 0.3903855
Авто- и кросс-корреляции нет
n = length(coef(ARMAmodelTANT))-1
x = var$p-n+1
rss = sum(ARMAmodelTANT$residuals[x:length(residuals(ARMAmodelTANT))]^2)
n1 = var$varresult$d1TATN$rank-1
ess = sum(var$varresult$d1TATN$residuals^2)
Fstat = ((rss-ess)/(n1-n))/(ess/(length(d1TATN)-2*n1-1))
Fstat
## [1] 4.565615
qf(0.99, df1=n1-n, df2=length(d1TATN)-2*n1-1)
## [1] 1.944196
pf(Fstat, n1-n, length(d1TATN)-2*n1-1, lower.tail=F)
## [1] 6.587105e-10
p-value = 6.587105e-10 => VAR(2) model превосходит ARMA model.
resvar = restrict(var, method = c("ser"), thresh = 0.8)
n1 = resvar$varresult$d1TATN$rank-1
if(n1<n+1){n1=n+1}
ess = sum(resvar$varresult$d1TATN$residuals^2)
Fstat = ((rss-ess)/(n1-n))/(ess/(length(d1TATN)-2*n1-1))
Fstat
## [1] 7.328169
qf(0.99, df1=n1-n, df2=length(d1TATN)-2*n1-1)
## [1] 2.257918
pf(Fstat, n1-n, length(d1TATN)-n-n1-1, lower.tail=F)
## [1] 2.29406e-12
Hosking(resvar, lags=1.5*var$p)
## lags statistic df p-value
## 15 25.81503 20 0.172008
LiMcLeod(resvar, lags=1.5*var$p)
## lags statistic df p-value
## 15 25.9219 20 0.1684068
F-критерий увеличился, при этом модель все еще лучше ARMA и не имеет авто- и кросс-корреляции остатков => превосходит модель без ограничений.
df = data.frame(d1TATN, d1SIBN, d1KMAZ)
VARselect(df, lag.max = 21, type = 'const')$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 10 2 1 10
var = VAR(df, p = 10, type = 'const')
Hosking(var, lags=1.3*var$p)
## lags statistic df p-value
## 13 26.19231 27 0.5079481
LiMcLeod(var, lags=1.3*var$p)
## lags statistic df p-value
## 13 26.48937 27 0.4915874
Авто- и кросс-корреляции нет
n = length(coef(ARMAmodelTANT))-1
x = var$p-n+1
rss = sum(ARMAmodelTANT$residuals[x:length(residuals(ARMAmodelTANT))]^2)
n1 = var$varresult$d1TATN$rank-1
ess = sum(var$varresult$d1TATN$residuals^2)
Fstat = ((rss-ess)/(n1-n))/(ess/(length(d1TATN)-2*n1-1))
Fstat
## [1] 3.456849
qf(0.99, df1=n1-n, df2=length(d1TATN)-2*n1-1)
## [1] 1.735509
pf(Fstat, n1-n, length(d1TATN)-2*n1-1, lower.tail=F)
## [1] 3.398337e-09
p-value = 3.398337e-09 => VAR(3) model превосходит ARMA model.
resvar = restrict(var, method = c("ser"), thresh = 0.7)
n1 = resvar$varresult$d1TATN$rank-1
if(n1<n+1){n1=n+1}
ess = sum(resvar$varresult$d1TATN$residuals^2)
Fstat = ((rss-ess)/(n1-n))/(ess/(length(d1TATN)-2*n1-1))
Fstat
## [1] 5.069553
qf(0.99, df1=n1-n, df2=length(d1TATN)-2*n1-1)
## [1] 1.915418
pf(Fstat, n1-n, length(d1TATN)-n-n1-1, lower.tail=F)
## [1] 5.914694e-12
Hosking(resvar, lags=1.3*var$p)
## lags statistic df p-value
## 13 34.46076 27 0.1530583
LiMcLeod(resvar, lags=1.3*var$p)
## lags statistic df p-value
## 13 34.73685 27 0.145658
F-критерий увеличился, при этом модель все еще лучше ARMA и не имеет авто- и кросс-корреляции остатков => превосходит модель без ограничений.
#модель с 2-мя переменными
df = data.frame(d1TATN, d1SIBN)
var = VAR(df, p = 10, type = 'const')
resvar = restrict(var, method = c("ser"), thresh = 0.8)
n = resvar$varresult$d1TATN$rank-1
x = 10-10+1
rss = sum(resvar$varresult$d1TATN$residuals[x:length(resvar$varresult$d1TATN$residuals)]^2)
#модель с 3-мя переменными
df = data.frame(d1TATN, d1SIBN, d1KMAZ)
var = VAR(df, p = 10, type = 'const')
resvar = restrict(var, method = c("ser"), thresh = 0.7)
n1 = resvar$varresult$d1TATN$rank-1
ess = sum(resvar$varresult$d1TATN$residuals^2)
Fstat = ((rss-ess)/(n1-n))/(ess/(length(d1TATN)-2*n1-1))
Fstat
## [1] 1.965555
qf(0.99, df1=n1-n, df2=length(d1TATN)-2*n1-1)
## [1] 2.52158
pf(Fstat, n1-n, length(d1TATN)-2*n1-1, lower.tail=F)
## [1] 0.04719159
F-критерий уменьшился, однако значение p-value на 5% уровне значимости говорит о том, что модель с тремя переменными корректна и лучше для объяснения d1TATN, чем с двумя.
df = data.frame(d1TATN, d1SIBN, d1KMAZ, d1NLMK)
VARselect(df, lag.max = 24, type = 'const')$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 8 1 1 8
var = VAR(df, p = 8, type = 'const')
Hosking(var, lags=1.2*var$p)
## lags statistic df p-value
## 9.6 29.88831 25.6 0.2543419
LiMcLeod(var, lags=1.2*var$p)
## lags statistic df p-value
## 9.6 30.20763 25.6 0.2416713
Авто- и кросс-корреляции нет
n = length(coef(ARMAmodelTANT))-1
x = var$p-n+1
rss = sum(ARMAmodelTANT$residuals[x:length(residuals(ARMAmodelTANT))]^2)
n1 = var$varresult$d1TATN$rank-1
ess = sum(var$varresult$d1TATN$residuals^2)
Fstat = ((rss-ess)/(n1-n))/(ess/(length(d1TATN)-2*n1-1))
Fstat
## [1] 3.148388
qf(0.99, df1=n1-n, df2=length(d1TATN)-2*n1-1)
## [1] 1.707842
pf(Fstat, n1-n, length(d1TATN)-2*n1-1, lower.tail=F)
## [1] 2.569228e-08
p-value = 2.569228e-08 => VAR(4) model превосходит ARMA model.
resvar = restrict(var, method = c("ser"), thresh = 0.6)
n1 = resvar$varresult$d1TATN$rank-1
if(n1<n+1){n1=n+1}
ess = sum(resvar$varresult$d1TATN$residuals^2)
Fstat = ((rss-ess)/(n1-n))/(ess/(length(d1TATN)-2*n1-1))
Fstat
## [1] 4.512438
qf(0.99, df1=n1-n, df2=length(d1TATN)-2*n1-1)
## [1] 1.864688
pf(Fstat, n1-n, length(d1TATN)-n-n1-1, lower.tail=F)
## [1] 5.348773e-11
Hosking(resvar, lags=1.2*var$p)
## lags statistic df p-value
## 9.6 34.52153 25.6 0.1118634
LiMcLeod(resvar, lags=1.2*var$p)
## lags statistic df p-value
## 9.6 34.8303 25.6 0.1052958
F-критерий увеличился, при этом модель все еще лучше ARMA и не имеет авто- и кросс-корреляции остатков => превосходит модель без ограничений.
#модель с 3-мя переменными
df = data.frame(d1TATN, d1SIBN, d1KMAZ)
var = VAR(df, p = 10, type = 'const')
resvar = restrict(var, method = c("ser"), thresh = 0.7)
n = resvar$varresult$d1TATN$rank-1
x = 10-8+1
rss = sum(resvar$varresult$d1TATN$residuals[x:length(resvar$varresult$d1TATN$residuals)]^2)
#модель с 4-мя переменными
df = data.frame(d1TATN, d1SIBN, d1KMAZ, d1NLMK)
var = VAR(df, p = 8, type = 'const')
resvar = restrict(var, method = c("ser"), thresh = 0.6)
n1 = resvar$varresult$d1TATN$rank-1
ess = sum(resvar$varresult$d1TATN$residuals^2)
Fstat = ((rss-ess)/(n1-n))/(ess/(length(d1TATN)-2*n1-1))
Fstat
## [1] -1.447424
qf(0.99, df1=n1-n, df2=length(d1TATN)-2*n1-1)
## [1] 4.617544
pf(Fstat, n1-n, length(d1TATN)-2*n1-1, lower.tail=F)
## [1] 1
Модель с 4-мя переменными хуже, чем модель с 3-мя переменными для объяснения d1TATN.
Конечная модель - VAR model с 3-мя переменными с ограничениями.
df_var = data.frame(d1TATN, d1SIBN, d1KMAZ)
var = VAR(df_var, p = 10, type = 'const')
resvar = restrict(var, method = c("ser"), thresh = 0.7)
Нас интересуют конкретные разрывы в 2020 году, а именно 31 января 2020 года - появились первые зараженные коронавирусом в России и 11 августа 2020 - регистрация в России первой в мире вакции от COVID-19. Поэтому рассмотрим разрывы на промежутке всего 2020 года.
Для новых “кусков” данных нам понадобятся новые модели ARMA для проверки структурных сдвигов. ln(250) = 5
Мы проводим Sup-F тест для поиска структурных разрывов, а также тест Bai Perron - для множественных разрывов.
d1TATN_sdvig = df_var$d1TATN[1261:1510]
ARMAmodelTATN_sdvig = Arima(d1TATN_sdvig, c(2,0,1), include.constant=TRUE, method = c("CSS-ML"))
Box.test(residuals(ARMAmodelTATN_sdvig), lag = 5, type = c("Ljung-Box"), fitdf = 3)
##
## Box-Ljung test
##
## data: residuals(ARMAmodelTATN_sdvig)
## X-squared = 0.10832, df = 2, p-value = 0.9473
Автокорреляции нет - модель годная.
d1TATN_sdvig = d1TATN_sdvig[1:length(d1TATN_sdvig)]
d1TATN_sdvig_l1 = c(0,d1TATN_sdvig[1:length(d1TATN_sdvig)-1])
d1TATN_sdvig_l2 = c(0,0,d1TATN_sdvig[2:length(d1TATN_sdvig)-2])
stat_TATN = Fstats(d1TATN_sdvig ~ d1TATN_sdvig_l1 + d1TATN_sdvig_l2, from = 0.1, to = NULL)
plot(stat_TATN, alpha = 0.01)
lines(breakpoints(stat_TATN))
a_TATN = breakpoints(stat_TATN)
a_TATN$breakpoints
## [1] 44
sctest(stat_TATN, type = "supF")
##
## supF test
##
## data: stat_TATN
## sup.F = 14.172, p-value = 0.05631
data$date[1261+44]
## [1] "2020-03-10"
Структурный разрыв есть на 10% уровне - 10 марта 2020 года.
d1SIBN_sdvig = df_var$d1SIBN[1261:1510]
ARMAmodelSIBN_sdvig = Arima(d1SIBN_sdvig, c(1,0,1), include.constant=TRUE, method = c("CSS-ML"))
Box.test(residuals(ARMAmodelSIBN_sdvig), lag = 5, type = c("Ljung-Box"), fitdf = 2)
##
## Box-Ljung test
##
## data: residuals(ARMAmodelSIBN_sdvig)
## X-squared = 2.3194, df = 3, p-value = 0.5088
Автокорреляции нет - модель годная.
d1SIBN_sdvig = d1SIBN_sdvig[1:length(d1SIBN_sdvig)]
d1SIBN_sdvig_l1 = c(0,d1SIBN_sdvig[1:length(d1SIBN_sdvig)-1])
stat_SIBN = Fstats(d1SIBN_sdvig ~ d1SIBN_sdvig_l1, from = 0.1, to = NULL)
plot(stat_SIBN, alpha = 0.01)
lines(breakpoints(stat_SIBN))
a_SIBN = breakpoints(stat_SIBN)
a_SIBN$breakpoints
## [1] 50
sctest(stat_SIBN, type = "supF")
##
## supF test
##
## data: stat_SIBN
## sup.F = 11.721, p-value = 0.05807
data$date[1261+50]
## [1] "2020-03-18"
Присутствует сдвиг на 10% уровне - 18 марта 2020 года.
d1KMAZ_sdvig = df_var$d1KMAZ[1261:1510]
ARMAmodelKMAZ_sdvig = Arima(d1KMAZ_sdvig, c(2,0,2), include.constant=TRUE, method = c("CSS-ML"))
Box.test(residuals(ARMAmodelKMAZ_sdvig), lag = 5, type = c("Ljung-Box"), fitdf = 4)
##
## Box-Ljung test
##
## data: residuals(ARMAmodelKMAZ_sdvig)
## X-squared = 2.6541, df = 1, p-value = 0.1033
Автокорреляции нет - модель годная.
d1KMAZ_sdvig = d1KMAZ_sdvig[1:length(d1KMAZ_sdvig)]
d1KMAZ_sdvig_l1 = c(0,d1KMAZ_sdvig[1:length(d1KMAZ_sdvig)-1])
d1KMAZ_sdvig_l2 = c(0,0,d1KMAZ_sdvig[2:length(d1KMAZ_sdvig)-2])
stat_KMAZ = Fstats(d1KMAZ_sdvig ~ d1KMAZ_sdvig_l1 + d1KMAZ_sdvig_l2, from = 0.1, to = NULL)
plot(stat_KMAZ, alpha = 0.01)
lines(breakpoints(stat_KMAZ))
a_KMAZ = breakpoints(stat_KMAZ)
a_KMAZ$breakpoints
## [1] 50
sctest(stat_KMAZ, type = "supF")
##
## supF test
##
## data: stat_KMAZ
## sup.F = 13.726, p-value = 0.06689
data$date[1261+50]
## [1] "2020-03-18"
Присутствует сдвиг на 10% уровне - 18 марта 2020 года.
На 10% уровне значимости наблюдаются стуктурные разрывы в марте (10 и 18 числа). Тогда начали вступать в силу первые серьезные ограничения из-за COVID-19, а заболеваемость в стране стремительно поднималась вверх.
Для сравнения моделей мы построили несколько разных прогнозов на год:
Прогноз на основе модели ARMA дает представление о колебании будущих цен на акции. Мы можем видеть, что цены везде идут вверх, но сильно медленнее поднимаются у акций Сургутнефтигаза и Аэрофлота.
forecast_TATN = forecast(ARMAmodelTANT, h=365)
plot(forecast_TATN)
forecast_SNGS = forecast(ARMAmodelSNGS, h=365)
plot(forecast_SNGS)
forecast_ALRS = forecast(ARMAmodelALRS, h=365)
plot(forecast_ALRS)
forecast_SIBN = forecast(ARMAmodelSIBN, h=365)
plot(forecast_SIBN)
forecast_NLMK = forecast(ARMAmodelNLMK, h=365)
plot(forecast_NLMK)
forecast_KMAZ = forecast(ARMAmodelKMAZ, h=365)
plot(forecast_KMAZ)
forecast_AFLT = forecast(ARMAmodelAFLT, h=365)
plot(forecast_AFLT)
Что касается прогноза на основе simple GARCH, он показывает прогноз именно изменения цен на акции. В целом, можно сказать, что прогнозы двух моделей дополняют друг друга.
prognoz_TATN = ugarchforecast(garch.TATN, n.ahead = 365)
sigmaforecast_TATN = c(garch.TATN@fit$sigma, sigma(prognoz_TATN))
plot(sigmaforecast_TATN, type = 'l')
prognoz_SNGS = ugarchforecast(garch.SNGS, n.ahead = 365)
sigmaforecast_SNGS = c(garch.SNGS@fit$sigma, sigma(prognoz_SNGS))
plot(sigmaforecast_SNGS, type = 'l')
prognoz_ALRS = ugarchforecast(garch.ALRS, n.ahead = 365)
sigmaforecast_ALRS = c(garch.ALRS@fit$sigma, sigma(prognoz_ALRS))
plot(sigmaforecast_ALRS, type = 'l')
prognoz_SIBN = ugarchforecast(garch.SIBN, n.ahead = 365)
sigmaforecast_SIBN = c(garch.SIBN@fit$sigma, sigma(prognoz_SIBN))
plot(sigmaforecast_SIBN, type = 'l')
prognoz_NLMK = ugarchforecast(garch.NLMK, n.ahead = 365)
sigmaforecast_NLMK = c(garch.NLMK@fit$sigma, sigma(prognoz_NLMK))
plot(sigmaforecast_NLMK, type = 'l')
prognoz_KMAZ = ugarchforecast(garch.KMAZ, n.ahead = 365)
sigmaforecast_KMAZ = c(garch.KMAZ@fit$sigma, sigma(prognoz_KMAZ))
plot(sigmaforecast_KMAZ, type = 'l')
prognoz_AFLT = ugarchforecast(garch.AFLT, n.ahead = 365)
sigmaforecast_AFLT = c(garch.AFLT@fit$sigma, sigma(prognoz_AFLT))
plot(sigmaforecast_AFLT, type = 'l')
Прогноз по VAR модели показывает прогноз целевой переменной – цен акций компании Татнефть - на основе предыдущих данных цен акций Газпром нефти и Камаз - влияющих на нее переменных.
На остатках:
df_var = data.frame(d1TATN, d1SIBN, d1KMAZ)
#VARselect(df_var, lag.max = 20, type = "const")
VAR_final = VAR(df_var, p = 10, type = "const")
prognoz_var = predict(VAR_final, n.ahead = 365)
plot(prognoz_var, nc = 1)
Прогноз основных данных:
Price_forARMA = c(data$TATN_p, forecast_TATN$mean)
Price_forARMA = Price_forARMA[1000:2130]
Price_upARMA = c(data$TATN_p, forecast_TATN$upper[,2])
Price_upARMA = Price_upARMA[1000:2130]
Price_lowARMA = c(data$TATN_p,forecast_TATN$lower[,2])
Price_lowARMA = Price_lowARMA[1000:2130]
plot(TATN_for, type = "l")
lines(TATN_up)
lines(TATN_low)
lines(Price_forARMA, col='red')
lines(Price_upARMA, col='red')
lines(Price_lowARMA, col='red')
Красная линия - прогноз по модели ARMA для цен на акции компании Татнефть.
Черная линия - прогноз на основе VAR(3) модели
На графике видно, что VAR прогноз немного ниже, чем арма, так как основано на днных трех переменных, а не одной, целевой.
Ожидаемая доходность: Ri = Rf + β(Rm−Rf)
Rf - безрисковая ставка доходности (доходность государственных облигаций - мы взяли годовые)
Rm - ожидаемая рыночная доходность (индекс Мосбиржи)
Rm − Rf - премия за инвестиционный риск
β - коэффициент чувтвительности (источник - Тинькофф Инвестиции)
Бета-коэффициенты: TATN_p = 1, SNGS_p = 0.99, ALRS = 0.93, SIBN = 0.86, NLMK = 0.27, KMAZ = 0.81, AFLT = 1.12
Мы считаем реальную доходность акций через процент изменения от дня ко дню, а доходность по модели CAPM – по формуле через линейную регрессию.
tatn = read_excel("stocks.xlsx", sheet = "TATN_p") %>% dplyr::select("Дата", "Изм. %") %>%
rename("date" = "Дата", "TATN_return" = "Изм. %")
sngs = read_excel("stocks.xlsx", sheet = "SNGS_p") %>% dplyr::select("Дата", "Изм. %") %>%
rename("date" = "Дата", "SNGS_return" = "Изм. %")
alrs = read_excel("stocks.xlsx", sheet = "ALRS") %>% dplyr::select("Дата", "Изм. %") %>%
rename("date" = "Дата", "ALRS_return" = "Изм. %")
sibn = read_excel("stocks.xlsx", sheet = "SIBN") %>% dplyr::select("Дата", "Изм. %") %>%
rename("date" = "Дата", "SIBN_return" = "Изм. %")
nlmk = read_excel("stocks.xlsx", sheet = "NLMK") %>% dplyr::select("Дата", "Изм. %") %>%
rename("date" = "Дата", "NLMK_return" = "Изм. %")
kmaz = read_excel("stocks.xlsx", sheet = "KMAZ") %>% dplyr::select("Дата", "Изм. %") %>%
rename("date" = "Дата", "KMAZ_return" = "Изм. %")
aflt = read_excel("stocks.xlsx", sheet = "AFLT") %>% dplyr::select("Дата", "Изм. %") %>%
rename("date" = "Дата", "AFLT_return" = "Изм. %")
imoex = read_excel("IMOEX.xlsx") %>% dplyr::select("Дата", "Изм. %") %>%
rename("date" = "Дата", "R_m" = "Изм. %")
imoex$date = ymd(imoex$date)
imoex$R_m = as.numeric(imoex$R_m)
bonds = read_excel("annual_bonds.xlsx") %>% dplyr::select("Дата", "Цена") %>%
rename("date" = "Дата", "R_f" = "Цена")
bonds$date = ymd(bonds$date)
capm = data %>% dplyr::select(date) %>% left_join(tatn, by = "date") %>% left_join(sngs, by = "date") %>%
left_join(alrs, by = "date") %>% left_join(sibn, by = "date") %>% left_join(nlmk, by = "date") %>%
left_join(kmaz, by = "date") %>%left_join(aflt, by = "date") %>%
left_join(bonds, by = "date") %>% left_join(imoex, by = "date")
capm$TATN_return = capm$TATN_return/100
capm$SNGS_return = capm$SNGS_return/100
capm$ALRS_return = capm$ALRS_return/100
capm$SIBN_return = capm$SIBN_return/100
capm$NLMK_return = capm$NLMK_return/100
capm$KMAZ_return = capm$KMAZ_return/100
capm$AFLT_return = capm$AFLT_return/100
capm$R_m = capm$R_m/100
capm$R_f = capm$R_f/100
capm = capm %>% mutate(risk_premium = R_m - R_f)
capm = capm %>% mutate(TATN_rp = 1 * risk_premium, SNGS_rp = 0.99 * risk_premium,
ALRS_rp = 0.93 * risk_premium, SIBN_rp = 0.86 * risk_premium,
NLMK_rp = 0.27 * risk_premium, KMAZ_rp = 0.81 * risk_premium,
AFLT_rp = 1.12 * risk_premium)
# считаем ожидаемую доходность
TATN_fit = lm(TATN_return ~ R_f + TATN_rp, data = capm)
SNGS_fit = lm(SNGS_return ~ R_f + SNGS_rp, data = capm)
ALRS_fit = lm(ALRS_return ~ R_f + ALRS_rp, data = capm)
SIBN_fit = lm(SIBN_return ~ R_f + SIBN_rp, data = capm)
NLMK_fit = lm(NLMK_return ~ R_f + NLMK_rp, data = capm)
KMAZ_fit = lm(KMAZ_return ~ R_f + KMAZ_rp, data = capm)
AFLT_fit = lm(AFLT_return ~ R_f + AFLT_rp, data = capm)
TATN_predict = predict.lm(TATN_fit)
SNGS_predict = predict.lm(SNGS_fit)
ALRS_predict = predict.lm(ALRS_fit)
SIBN_predict = predict.lm(SIBN_fit)
NLMK_predict = predict.lm(NLMK_fit)
KMAZ_predict = predict.lm(KMAZ_fit)
AFLT_predict = predict.lm(AFLT_fit)
capm$TATN_predict = TATN_predict
capm$SNGS_predict = SNGS_predict
capm$ALRS_predict = ALRS_predict
capm$SIBN_predict = SIBN_predict
capm$NLMK_predict = NLMK_predict
capm$KMAZ_predict = KMAZ_predict
capm$AFLT_predict = AFLT_predict
Красная линия - реальная доходность Синия линия - доходность по модели CAPM
date_d = c(1:1765)
capm = capm %>% mutate(date_d = date_d) %>% filter(date_d > 1399)
capm$date = ymd(capm$date)
ggplot(capm) +
geom_line(aes(x = date, y = TATN_return), color = "red") +
geom_line(aes(x = date, y = TATN_predict), color = "blue") +
ggtitle('График доходности по модели CAPM и реальной доходности \n Татнефть') +
xlab('Дата') +
ylab('Доходность') +
theme_bw()
ggplot(capm) +
geom_line(aes(x = date, y = SNGS_return), color = "red") +
geom_line(aes(x = date, y = SNGS_predict), color = "blue") +
ggtitle('График доходности по модели CAPM и реальной доходности \n Сургутнефтегаз') +
xlab('Дата') +
ylab('Доходность') +
theme_bw()
ggplot(capm) +
geom_line(aes(x = date, y = ALRS_return), color = "red") +
geom_line(aes(x = date, y = ALRS_predict), color = "blue") +
ggtitle('График доходности по модели CAPM и реальной доходности \n Алроса') +
xlab('Дата') +
ylab('Доходность') +
theme_bw()
ggplot(capm) +
geom_line(aes(x = date, y = NLMK_return), color = "red") +
geom_line(aes(x = date, y = NLMK_predict), color = "blue") +
ggtitle('График доходности по модели CAPM и реальной доходности \n Газпром нефть') +
xlab('Дата') +
ylab('Доходность') +
theme_bw()
ggplot(capm) +
geom_line(aes(x = date, y = KMAZ_return), color = "red") +
geom_line(aes(x = date, y = KMAZ_predict), color = "blue") +
ggtitle('График доходности по модели CAPM и реальной доходности \n Новолипецкий металлургический завод') +
xlab('Дата') +
ylab('Доходность') +
theme_bw()
ggplot(capm) +
geom_line(aes(x = date, y = AFLT_return), color = "red") +
geom_line(aes(x = date, y = AFLT_predict), color = "blue") +
ggtitle('График доходности по модели CAPM и реальной доходности \n Камаз') +
xlab('Дата') +
ylab('Доходность') +
theme_bw()
ggplot(capm) +
geom_line(aes(x = date, y = SIBN_return), color = "red") +
geom_line(aes(x = date, y = SIBN_predict), color = "blue") +
ggtitle('График доходности по модели CAPM и реальной доходности \n Аэрофлот') +
xlab('Дата') +
ylab('Доходность') +
theme_bw()
Таким образом, мы построили графики сравнения реальной доходности и доходности, прогнозируемой моделью. В целом, можно увидеть, что модель довольно четко описывает реальную доходность, имея лишь меньше выбросов.